Modelling charge self-trapping in wide-gap dielectrics: 
Localization problem in local density functionals 
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We discuss the adiabatic self-trapping of small polarons within the density functional theory 
(DFT). In particular, we carried out plane- wave pseudo-potential calculations of the triplet exciton 
in NaCl and found no energy minimum corresponding to the self-trapped exciton (STE) contrary 
to the experimental evidence and previous calculations. To explore the origin of this problem we 
modelled the self-trapped hole in NaCl using hybrid density functionals and an embedded cluster 
method. Calculations show that the stability of the self-trapped state of the hole drastically depends 
on the amount of the exact exchange in the density functional: at less than 30% of the Hartree-Fock 
exchange, only delocalized hole is stable, at 50% - both delocalized and self-trapped states are stable, 
while further increase of exact exchange results in only the self-trapped state being stable. We argue 
that the main contributions to the self-trapping energy such as the kinetic energy of the localizing 
charge, the chemical bond formation of the di-halogen quasi molecule, and the lattice polarization, 
are represented incorrectly within the Kohn-Sham (KS) based approaches. 
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I. INTRODUCTION 

Despite considerable progress in studies of self-trapped 
excitons and polarons, the dynamics of early stages of 
self-trapping in specific systems is still poorly under- 
stood. The conceptual difficulty primarily lies with the 
fact that the quasi particle in question undergoes a tran- 
sition between the free (delocalized) state and the lo- 
calized one. In case of small polarons this transition is 
also associated with a substantial local lattice distortion. 
Within an adiabatic approximation one can describe self- 
trapping in terms of a potential energy surface (PES) 
connecting free and localized states. Then a microscopic 
model of self-trapping will involve a characterization of 
this PES, i.e. relevant atomic coordinates, relative ener- 
gies of the free and self-trapped states, the energy barrier 
(if any) between them, as well as spectroscopic properties 
of free and self-trapped species. 

Significant progress has been achieved in understand- 
ing of the conditions for polaron and exciton localiaajtk 
and self-trapping, as summarized in recent reviewsatrL 
It has been suggested by RashbaO that in three dimen- 
sional dielectric crystals, the free and self-trapped forms 
of excitons may coexist, thus implying an energy bar- 
rier separating the two states. The height of this bar- 
rier affects the dynamics and characteristic time of self- 
trapping process. However, atomistic modelling has only 
been successful in calculating the structure and proper- 
ties of strongly localized systems and never in modelling 
transitions between delocalized and localized states. One 
should note that the barrier for self-trapping, if it exists, 
is not likely to exceed a few tenths of an electron voltfl. 
Therefore, its experimental verification and theoretical 
calculation is extremely challenging. 

It has been anticipated that further development 
of quantum mechanical techniques, especially Density 
Functional Theory (DFT), will allow one to close this 
gap and achieve predictive modelling of self-trapping or 
defect induced trapping process. However, recent at- 



tempts to calculate even well established models of small 
radius polarons have failed unexpectedly. Calculations 
of the triplet self-trapped exciton in NaCl using plane- 
wave DFT in the Generalized Gradient Approximation 
(GGA) predicted no stable statetl in direct contradiction 
with the experimental cvidenceo and previous Hartree- 
Fogk calculationsD'El. On the other hand, Perebeinos et 
al.u using plane- wave local spin density (LSDA) approach 
predicted the existence of a marginally stable STE in 
this system, albeit highej-pip energy than the free exci- 
ton state. Song et a/HUtS applied plane wave DFT to 
the triplet exciton in a-quartz and similarly found the 
free exciton state to be more stable than the localized 
one. The delocalized solution has also been found in 
LSDA calculations of the hole trapping at the Li _center 
in MgCO. Pacchioni et aZo and Laegsgaard et aZEj con- 
sidered the hole trapped at an Al impurity in silica. Both 
groups concluded that the predicted structure of this de- 
fect strongly depends on the density functional used in 
the calculations: local and GGA functionals predicted 
only the delocalized hole to be stable, again in contra- 
diction with the experiment. Importantly, the non-local 
density functionals predict more localized states for this 
system. 

An apparent bias of DFT calculations towards the 
delocalized electronic states was attributed to the self- 
interaction error inherent in the local or semilocal GGA 
type approximations, which are central to the Kohn- 
Sham methodllS. It is unclear, however, to what extent 
the local approximation affects qualitative results (local- 
ized versus delocalized states), and what is the role of 
other factors in the calculation, such as boundary con- 
ditions, basis set completeness, and pseudopotential ap- 
proximation. 

In this paper, we consider a triplet exciton and a hole 
in the archetypal ionic insulator NaCl. These defects 
have been studied extensively both experimentally and 
theoreticallyoiJ : and thus provide good test systems. 
Previous calculations were carried out mainly in small 



embedded cluster models using the Hartree-Fock method 
and therefore were unable to treat delocalized states and 
take full account of the electron correlation. We would 
like to model a transition from delocalized to self-trapped 
exciton state and for this purpose use the plane-wave 
DFT method. We analyse the effect of size of a periodic 
supercell and the related question of spurious multipole 
interactions in this system and conclude that no stable 
self-trapped state for the exciton or a hole is predicted 
within the GGA DFT framework, once these factors are 
eliminated. To separate the localization problem from 
the effects of periodic boundary conditions (PBC) we 
then consider a hole within the embedded cluster ap- 
proach and find that the non-local contribution to the 
exchange interaction is decisive in the description of self- 
trapped states. 

The paper is organized as follows. In the next Section 
we give a brief account of the microscopic models of the 
self-trapped exciton and hole in alkali halides. Next we 
outline the details of the periodic plane- wave and embed- 
ded cluster procedures used. The results of calculations 
for an exciton and a hole are presented respectively in 
Sections IV and V, followed by discussion in section VI. 







II. SELF-TRAPPED EXCITON AND HOLE IN 
ALKALI HALIDES: BACKGROUND 

Most of the alkali halides at normal pressure and tem- 
perature assume the face centered cubic structure (Fig. 
[j]a) (except cesium halides which are simple cubic). Self- 
trapped excitons (STE) and holes in these crystals ex- 
hibit the features of both molecular and dielectric polaron 
character. In particular, the self-trapped hole (otherwise 
called a 14-center) is known to be localized on two adja- 
cent halogen ions forming a X^ molecular ion (Fig. [j]d). 
There is no experimental evidence pointing to existence 
of a barrier for the hole self-trapping in alkali halidesBEEl, 
suggesting, that the free hole is unstable. In particular, 
LushchikEj et al. estimated that the mean free path of 
a free hole in NaCl before self-trapping does not exceed 
30 a,Q (ao is the shortest Na-Cl distance in the perfect 
lattice). The hole self-trapping energy, which is the dif- 
ference between the energies of the fully delocalized and 
localized states, is larger than the activation energy for 
the Vfe center diffusion (~ 0.4 eV in NaCl)E3. Therefore, 
it is expected that the description of such a deep state is 
well within the reach of the DFT theory. 

The model and the geometric structure of the 14-center 
in alkali halides has been first proposed on thc-hasis of the 
analysis of the electron spin resonance dataE2lEi:-ax 
refined in numerous theoretical calculationstjt 
The structure and stability of the 14-center is primar- 
ily determined by the chemical bond formation between 
the two halogen ions, assisted by the lattice polarization. 

Self-trapped triplet excitons in alkali halides can be 
formed either directly by excitation in the exciton band 
or as a result of recombination of electron with a self- 



FIG. 1: Schematic representations of the perfect rock-salt 
structure (a); of the triplet self-trapped exciton (b), the clos- 
est separated F-H pair (c), and the 14-center (d) in an alkali 
halide crystal. The marked atoms were constrained during 
the optimization procedure. 



trapped hole, 14-center. It is currently accepted that the 
STE in alkali halides consists of a hole localized on the 
X% molecular ion, and an electron localized in its vicin- 
ity. The so-called on-center and off-center configurations 
of triplet STE are considered in the literatureo. In the 
on-center configuration the X^ molecular ion equally oc- 
cupies two nearest anion sites (as in Fig. 0d), with the 
electron cloud symmetrically localized around it and has 
the D2h point symmetry group. It is therefore called a 
(Vk + e) configuration. Williams and Songc3 suggested 
that in some alkali halides at least, the on-center config- 
uration of the STE is unstable relative to an off-center 
displacement of the X^ ion along the < 110 > crystal- 
lographic direction towards one of the anion sites (Fig. 
|l|b). The center of mass of the electron component in 
this configuration is localized near the other (vacant) an- 
ion site, thus giving the Civ symmetry to the STE. This 
model of the triplet STE in alkali halides has been further 
developed using effective potential and embedded cluster 
ab initio Hartree-Fock calculations, as reviewed by Song 
and WilliamsB. The STE in NaCl has also been studied 
using ab initio embedded cluster methods by PuchiroEII 
et at, who confirmed the existence of the off-center STE 
configuration and demonstrated the importance of the 
electron correlation in determining the properties of the 
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STE. 

Triplet STE created by irradiation may subsequently 
annihilate (radiatively or non-radiatively) restoring the 
perfect lattice, or decay into a Frenkel pair of lattice de- 
fects: an electron in an anion vacancy (F-center) and an 
interstitial halogen atom in the form of an AT^ molecular 
ion occupying one halogen site, the H-center, (Fig. |l|c). 
The energy barrier between the off-center STE and the 
closest F-H pair has been calculated in alkali chlorides 
and bromides by Song et|-al.E3 and also recently in KBr 
by Shluger and Tanimuraa. The dissociation barrier for 
the STE is found-to be very low (0.07 eV for NaCl and 
0.03 eV for KBr)N. The activation .eaergy for the STE 
diffusion in NaCl is somewhat largeinj ~ 0.13 eV. 

At the outset of this work we hoped to use the plane 
wave periodic DFT methods to treat both delocalized 
and localized states and to study the adiabatic potential 
for self-trapping of triplet exciton in NaCl. The results 
presented below highlight the limited applicability of the 
existing functionals to the localization problem. 



III. CALCULATION PROCEDURE 

A. Periodic DFT calculations 

In order to model the triplet exciton self-trapping we 
consider the adiabatic coordinate corresponding to a dis- 
placement, A, of one chlorine ion (hereafter designated 
as CI*) along the < lUL > crystalline axis with all 
other surrounding atomstJ relaxed at each position of 
the CI*. The displacement A is measured in units of 
1/8 of the plane diagonal (ao\/2). By virtue of construc- 
tion, the specified adiabatic coordinate connects the free 



la) with the STE config- 
lb), and the closest F-H 



exciton state (A = 0) (Fig. 
uration (A ~ 4 - 5) (Fig. 
pair separated by one lattice anion, A ~ 8) (Fig. |l|c). 
The considered potential energy surface (PES) allows for 
both on- and off-center forms of STE as well as for the 
transient .one-centre excitons proposed by Shluger and 
Tarumuran. However, it should be emphasised, that an 
existence of the potential energy barrier along this coor- 
dinate does not rule out a possibility of alternative self- 
trapping paths with still lower energy barriers. 

The PES along the coordinate A was calculated us- 
ing the plane-wave density, functional approach imple- 
mented in the VASP codeEj with the GGA density func- 
tional due to Ppdew and WangEil and the ultra-soft 
pseudopotentialsEa as supplied by Kresse and Hafheio. 
The self-consistent triplet excited state was modelled by 
constraining the total spin of the system to S=l, so as to 
obtain the lowest energy state of given multiplicity. 

To study the effect of the super-cell size, the exciton 
PES was calculated in unit cells containing 32, 64, 108 
and 144 ions. The details of the super-cells are given 
in Table | together with the sets of k-points used and 
number of atoms explicitly relaxed. At each position of 
CI* the specified ions were relaxed, so that the maximum 



TABLE I: The supercells used in the exciton calculations, 
specified by the number of atoms, Cartesian translation vec- 
tors (given in units of the shortest Na-Cl distance a ) and the 
number of irreducible k-points used in the calculations. The 
last column shows the value of the optimized lattice constant 
do for the ground and excited states. 



Atoms ■- 


| Translation vectors (ao) 


No k-points 


ao 


A 


(relaxed)E 






s=0 


s=l 


2 


(1 1 0), (1 1), (0 1 1) 


35 


2.845 




32 (26) 


(4 0), (0 4 0), (2 2 2) 


4 


2.83 


2.87 


64 (48) 


(4 0), (0 4 0), (0 4) 


4 


2.84 


2.85 


108 (28) 


(3 3 0), (3 3), (0 6 6) 


3 


2.83 


2.84 


144 (46) 


(6 0), (0 6 0), (0 4) 


1 


2.84 


2.84 



force acting on any individual atom did not exceed 0.04 
eVf A cB. Most of the calculations were performed with 
an energy cut-off of 285 eV. 

As shall be seen below, the calculated PES for the 
exciton self-trapping does not display a minimum for the 
STE. Furthermore, the shape of PES strongly varies with 
the size of the supercell. In order to separate the effects 
of the boundary conditions and density functionals, we 
considered the formation of the V^-center in NaCl using 
an embedded cluster method implemented in the GUESS 
code 



B. Embedded cluster calculations 

In this method, a crystal is represented by a large finite 
nano-cluster, that is divided into several regions. Spher- 
ical region I at the centre of the nano-cluster includes a 
quantum-mechanically treated cluster (QC) surrounded 
by interface ions and a region of classical shell model 
ionscS. The remaining part of the nano-cluster is rep- 
resented by classical non-polarizable ions. All quantum 
mechanical, interface and classical ions (both cores and 
shells) in region I are allowed to relax simultaneously in 
the course of geometry optimization. Ions outside re- 
gion I remain fixed and provide an accurate electrostatic 
potential within region I. This approach allows one to 
take into account the defect-induced lattice polarization 
of a crystal region containing a few hundred atoms. An 
account of lattice polarisation outside region I can be ex- 
tended to infinity using a polarizsble continuum model 
and the Mott-Littleton approacbEl In this approxima- 
tion the polarization energy is proportional to the square 
of the excess charge in the lattice. In the case of the 
Vfe-center, the charge remains constant, so the Mott- 
Littleton contribution cancels out in any energy differ- 
ences. Therefore this term was neglected in the current 
calculations. 

The.[d£scribed scheme is implemented in the GUESS 
codeoEj, which provides an interface between the 
GAUSSIAN packageSa for calculation the electronic 
structure of the QC and a code providing the shell model 
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FIG. 2: Illustration of the embedded cluster model used in the 
Vfc-center calculation (only region 1 is shown). The quantum 
cluster is shown enlarged. The small white balls represent the 
sodium interface atoms 



representation of the rest of the crystal. The total energy 
and the electronic structure of the QC is calculated by 
solving standard Hartree-Fock or Kohn-Sham equations 
which include the matrix elements of the electrostatic po- 
tential due to all classical point charges in regions I and II 
computed on the basis functions of the cluster. Further 
details for the total-energy evaluation in this approach 
are given elsewherea. 

A cubic nano-cluster containing 8000 (20x20x20) ions 
was used in these calculations. Spherical region I of ra- 
dius 11.5 A included 300 ions and a stoichiometric QC 
comprising 22 ions (Fig. ||). 

In our scheme, the charge density in QC interacts with 
the classical atoms only electrostatically, which results in 
artificiai-polarization by the classical cations nearest to 
the QCE3. To avoid this, an interface region between the 
QC and classical ions in region I was introduced. In the 
present model it includes 20.JjJa ions represented by ef- 
fective core pseudopotentialsc2l, and a single s-type basis 
function (see Fig. ^). Therefore all QC anions are fully 
coordinated by either quantum or interface ions. For 
the quantum cluster wa.have used an all electron Gaus- 
sian type atomic basisCJ optimized for the NaCl crystal 
(8 - 511G for Na and 86 - 311G for CI). 

The shell model ions in region I interact via pair po- 
tentials due to Catlow et aLc3, which were slightly mod- 
ified to represent more accurately the static dielectric 
constant. Modified potentials give a static equilibrium 
Na — CI distance of 2. 79 A and fairly well reproduce the 
elastic, dielectric and phonon properties of the perfect 
NaCl crystal. 

The embedded cluster calculations were carried out 
using three different density functionals: the Hartree- 
Fock, Becke hybrid three parameter (B3), and so-called 
half-and-half functionals (BH&H)[3. The latter two in- 
corporate 20% and 50% of the Hartree-Fock exchange, 
respectivelyE^. The hybrid exchange functionals were 
used in conjunction with the Lee, Yang and Parr cor- 



relation functional (LYP)Q 

Consistency of the interactions between the classical 
and quantum regions was tested by calculating an ideal 
QC embedded into an ideal classical lattice. The opti- 
mized Na — CI distances inside the QC and those at the 
cluster boundary differ from the Na — Cl distances in the 
rest of the nano-cluster by ca. 1% and 3.5% respectively. 



IV. EXCITON SELF-TRAPPING IN NACL 



Exciton self-trapping is associated with a strong lattice 

:ation is 
1451 To ac- 



relaxation. It has been suggested that this re 
associated with the significant volume change 
count for this effect, we optimized the lattice constant 
of the supercells in both the ground and excited triplet 
states using the VASP code. The results are presented 
in (Table |). The calculated ground state equilibrium 
lattice constant ao — 2.84 A is ~ 2% larger than the ex- 
perimental value, extrapolated to K, of 2.79 A . We also 
observe a slight increase of the perfect lattice constant in 
the triplet excited state in smaller supercells. 

Next, we calculate the potential energy surface for the 
triplet exciton by displacing a selected CI* ion along the 
< 110 > crystallographic axis and relaxing all the sur- 
rounding atoms at each CI* position (Fig. |l|a,b,c). The 
PES for the triplet state calculated in different super- 
cells is shown in Figure ^. Additional information on 
the lattice relaxation and spin density changes along the 
CI* displacement coordinate can be found on our web 
pagec3. It is seen that the obtained PES varies signifi- 
cantly with the size and shape of the supercell, in which 
it was calculated. In all calculations the free exciton con- 
figuration is predicted to be the lowest energy state, as 
reported previouslyBB. The off-center STE configuration 
corresponds approximately to (A ~ 5). This state is pre- 
dicted to be marginally stable only in the smallest suppjr- 
cell used (identical to the one used by Perebeinos et alB). 
We also observe that the energy difference between the 
free exciton and the expected self-trapped state increases 
with the size of the supercell. 

We have singled out most probable causes of the signif- 
icant dependence of the calculated PES on the size and 
shape of different supercells: i) defect volume change; ii) 
spurious electrostatic interaction between the supercells; 
iii) basis set completeness, and iv) pseudopotential ap- 
proximation. In the rest of this Section we discuss these 
aspects. 

To mimic the effect of the volume change associated 
with the self-trapping, we optimized the volume of the 
144 ion supercell at each point of the adiabatic surface. 
The equilibrium lattice constant, ao has increased only 
slightly: from 2.840 A for a perfect lattice (A = 0) to 
2.848 A at the expected STE state (A = 5). The total 
energy of the STE has lowered by ~ 0.1 eV and the 
energy minimum for the STE did not exist. 

As can be seen in (Fig. |l|b), the charge distribution in 
off-center STE is associated with a substantial dipole mo- 
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FIG. 3: The potential energy surface for the triplet exciton as 
described in the text. The energies are plotted with respect 
to the energies of the triplet state of a perfect lattice (zero 
coordinate). A ~ 5 corresponds to the expected minimum for 
the STE (Fig. 0b), A = 8 corresponds to the next nearest- 
neighbour F-H pair (Fig. |l|c). The dipole corrected (solid 
line) and uncorrected (dashed line) energies are shown for the 
supercell of 144 ions. 

merit. The strong dependence of the PES on the supercell 
size suggests that a significant contribution to the total 
energy can arise from the interaction between periodic 
dipoles persisting even for comparatively large supercells. 
The dipole-dipole interaction energy resulting from the 
super-lattie« of dipoles can be evaluated according to the 
expressionED: 

2G 3 1 

Edip = 7T1=P + 7;PaT a pPf3, (1) 
oyTT Z 

where Einstein notation is adopted for the sums over re- 
peated indices, Pi denotes the Cartesian component of 
the dipole vector P, G is the Ewald parameter, and 

H ° z^O 

Here the first sum is over the reciprocal lattice vectors 
R g and the second - over the direct lattice vectors Ri 
with the zero vector (i = 0) being excluded to avoid 
self-interaction. In addition, H a p(y) = —6 a pH{y) + 
y a yp/y 2 m{y)+4/^e.y 2 } and H(y) = 2^y- 2 e^ + 
erfc(y) /y 3 . Equation (|l|) is a generalization of the Makov 
and Payne formulae^ for the non-cubic supercells. 

To quantify the dipole-dipole interaction energy the 
dipole moment of the supercell needs to be calculated. 
For a neutral periodic system this, must be defined as an 
invariant of a coordinate systcmc3: 

P(r) = n- 1 f rp(r)rf 3 r + fi" 1 f r [nll(r)] ds, (3) 

J cell J surface 



where the first integral is taken over the volume of the su- 
percell, and the second over the supercell surface, fl is the 
volume of the supercell, p(r) is a charge density (includ- 
ing both electronic and atomic cores), n is an outward 
surface-normal unit vector, and II is a local polarization 
defined via the equation: 

Vn(r) = -p(r). (4) 

It should be noted, that within the PBC only the whole 
sum in Eq. ^| is an invariant with respect to the co- 
ordinate origin. However, the second integral there can 
be made small by an appropriate choice of a coordinate 
system. Then and only then the dipole moment of the 
supercell can be correctly estimated by the first integral 
in Eq. (^). To minimize the contribution of the surface 
integral we chose the supercell as a Wigner-Seitz cell, 
where the origin of the charge density grid is located not 
on an atomic site, but on a volume interstitial position 
in the face centered cubic lattice. 

The Cartesian components of the calculated dipole mo- 
ment for each atomic configuration of the exciton along 
the adiabatic coordinate, A, in the supercell of 144 atoms 
are shown in Fig. ^. The dipole moment varies non- 
monotonically along the A coordinate reaching a max- 
imum near the expected STE configuration. Thus, the 
PES corrected for the dipole-dipole interaction may in 
principle contain a local minimum for the STE. To ex- 
plore this possibility we calculated the value of the dipole 
energy (Eq. |l|) in different unit cells for the expected 
STE configuration (A = 5). For that we assumed that 
calculations in different unit cells yield approximately the 
same charge density for the STE and used the value for 
the dipole moment, P, calculated in the 144 supercell for 
A = 5 (Fig. ||). The obtained values of Edi P are equal 
to -4.85, -2.64, 0.18, and -0.27 eV for the 32, 64, 108, 
and 144 ion supercells, respectively, and are negative ex- 
cept for the 108 ion supercell. This difference is caused 
by the fact that this supercell has a rhombohcdral shape 
whereas the other cells are rectangular. 

From the tendency shown in Fig. ^, we expect that the 
Edi P (A = 5) values given above represent the largest cor- 
rections for different supercells. When subtracted from 
the total energy, E(A) — E < n p (A), they in most cases en- 
hance even further the difference between the free and 
off-center STE energies. The full PES corrected by the 
dipole-dipole contribution, E(A) — Edi p (A), was calcu- 
lated only for the largest 144 ion supercell where the 
correction is the smallest (see Fig. ||). Although the 
shape of the PES changes slightly, it does not show a 
minimum for any displacement A. It is also clear that 
no energy minimum for the STE will occur in the 32 ion 
supercell after the dipole correction is taken into account. 
Our calculations clearly demonstrate that spurious elec- 
trostatic interaction in small unit cells may result in an 
artificial energy minimum for the STE. We suggest that 
the marginally stable STE configuration in the 32 atoms 
supercell reported by PerebeinosB et al. and also found 
in our calculation is precisely of this origin. 
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FIG. 4: The electric dipole moment as a function of the adi- 
abatic coordinate A, remaining details as in Fig. ^ 

We should note that the Edi P calculated using Eq. (jl]) 
overestimates the electrostatic interaction between the 
supercells in PBC due to various approximations in Eq. 

® : n 

1. The point dipole approximation assumed in Eq. (jfj), 
breaks down whenever the dipole length is comparable 
with the supercell size. (Hence unrealistically large dipole 
energies are obtained in small supercells). 

2. Eq. (0) does not account for a dielectric screening. 

3. The higher multipole contributions (dipole- 
quadrupole etc.) would generally have the sign oppo- 
site to E,n p and may reduce the overall effect of spurious 
multipole interactions. 

4. Due to a strong multipole interaction in smaller 
supercells, the charge distributions, and therefore the 
dipole moments, predicted in different supercells may dif- 
fer significantly. 

In spite of these remarks we believe that further re- 
finement of the calculations will not affect the main con- 
clusion - no stable state for the STE is predicted by the 
GGA DFT approach. 

The reasons for the absence of a STE energy minimum 
in the local DFT approximation are not obvious. Several 
factors could result in the bias towards the delocalized 
solution. These include, for example, the choice of en- 
ergy cut-off and k-point sampling. To check these issues, 
we compared calculations with cut-off energies of 220 cV 
and 285 eV and observed no significant difference in the 
STE relaxation. The k-point sampling is likely to affect 
more the energies of delocalized states. In fact we do ob- 
serve a slight increase of the free exciton energy with the 
number of irreducible k-points. This reduces the energy 
difference between the free and the self-trapped exciton, 
though the minimum for the STE still does not arise. An- 
other uncontrolled factor in the calculation is the pseu- 
dopotential approximation for core electrons. In particu- 



lar, we employed a large core (ls 2 2s 2 p 6 ) pseudopotential 
for sodium atoms and the polarization of cations was vir- 
tually neglected in our calculations. To examine the role 
of cation polarization, we carried out the STE calcula- 
tions with the Helium core pseudopotential for Na atoms 
for selected configurations in the 64 atoms supercell. This 
also did not affect our qualitative conclusions - no energy 
minimum was obtained for the STE. 

Our analysis suggests that the problem of localization 
lies not with the technical details of the calculations but 
rather more fundamentally with the approximations in- 
volved in the GGA density functional. In particular, es- 
sentially uncontrolled self-interaction error in the local 
density approximation is the most likely cause. How- 
ever, a detailed characterisation of the involved energy 
terms is difficult for the exciton, which comprises two 
strongly interacting spins with the energies near the top 
of the valence band and the bottom of the conduction 
band respectively. It must be emphasised at this point 
that holes in alkali halides do self-trap in the form of 
T4-centres, while the electrons do not. Hence, the hole 
self-trapping is also fundamental in the exciton localiza- 
tion. At the same time, self-trapping of the hole is a 
simpler problem in that it involves localization of a sin- 
gle spin. However, being charged, a hole presents some 
difficulty for studying within the PBC. In particular, we 
obtained no self-trapping of the hole within plane wave 
DFT. Therefore, we resorted to modelling the hole within 
the embedded cluster method, that has the following ad- 
vantages: I. The atomic- type basis set ensures no bias 
towards the delocalized states; 2. The lattice polarisation 
effects paramount for charged polar systems are fully ac- 
counted for; 3. The use of the atomic basis sets allows 
the straightforward application of a number of local and 
non-local density functionals, so that the properties of 
various functionals can be studied. 



V. HOLE TRAPPING IN NACL 

The calculations of the hole (a total charge of the sys- 
tem is equal to +1) were performed within the embedded 
cluster model discussed in Section IIIB (Fig. 2). We cal- 
culated the total energy of the system as a function of the 
distance between the two CI ions equally displaced along 
the < 1 10 > axis towards each other from their respec- 
tive perfect lattice sites (see Fig. |d). All other ions m 
region I were relaxed simultaneously to their equilibrium 
positions for each Cl-Cl separation. 

The potential energy surfaces calculated with differ- 
ent density functionals are depicted in Fig. 0. The zero 
energy in the top panel corresponds to the singly ion- 
ized perfect lattice calculated within each method. We 
observe that the shape of PES qualitatively depends on 
the amount of the Hartree-Fock exchange included in the 
density functional. The B3LYP functional predicts one 
global minimum corresponding to a free hole, the BH&H 
functional predicts two minima separated by a marginal 
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Cl 2 distance (A) 

FIG. 5: The adiabatic surface for the 14-center (top panel), 
calculated with various density functionals. At each given 
CI 2 distance, all other atoms in the system were relaxed to 
the minimum energy geometry. Zero energy corresponds to 
the perfect crystal configuration (the extended hole). The 
middle panel shows the dependence of the Mulliken spin pop- 
ulation on the CI2 molecule as a function of the same adia- 
batic coordinate. The bottom panel shows the relaxation of 
the pair of Na atoms adjacent to the Vfc- center. 



barrier, while the Hartree-Fock method gives a single 
energy minimum corresponding to the self-trapped hole 
(V^-center) (top panel of Fig. ||). Furthermore, the en- 
ergy difference between the free and self-trapped states 
ranges between -0.2 eV for B3LYP and 1.5 eV for the 
HF. 

Different functionals also predict distinctly different 
spin localization. The analysis of the Mulliken spin pop- 
ulation on the CI2 molecule (middle panel in Fig. ||) re- 
veals that the extent of spin localization increases faster 
with decreasing Cl^ distance when the amount of exact 
exchange is larger. In particular, the Hartree-Fock cal- 
culations predict the spin almost fully localized on Cl% 
at a mere 0.05 A displacement of two CI ions from their 
respective lattice sites. The B3LYP functional, on the 
contrary, predicts very gradual spin localization with the 
decrease of the Cl-Cl distance. The variance in localiza- 
tion of spin in different functionals also affects the lattice 
relaxation. The bottom panel in Figure [5] shows the dis- 
tance between the two relaxed Na ions adjacent to the 
CI2 molecule as a function of the Cl^ distance. It is seen 
that weak spin localization predicted in the B3LYP func- 
tional causes much smaller cation relaxation than that 
obtained in the Hartree-Fock method. 



VI. DISCUSSION 

Our results highlight severe problems with applying 
density functional theory to the self-trapping problem. 
The LSDA and GGA density functionals yield solutions 
biased towards the delocalized states to such an extent 
that no stable self-trapped state is predicted for either a 
hole or an exciton, contrary to experimental evidence. 

Let us emphasise that the degree of localisation of a 
polaron is closely linked to the properties of the single 
particle spectrum for the problem. For instance, spin 
density of the T4-center is represented almost entirely by 
the single-particle density of the last occupied orbital in 
a majority spin, i[; s (r'). The degree of localization of 
ip s (r') is related to its energy, e s , or more precisely, to 
the energy split of the e s from the corresponding band 
edge. In the KS method, e s is defined asliS 

( s =< ip s \ - ^V 2 +« e //(r)|V» >, (5) 

where 

v eff (r) = v(r) + J y^dr' + v xc (r). (6) 

Here, v(r) represents the external potential of the atomic 
cores, the second term is the classical electrostatic poten- 
tial due to electrons of density p(r), (including the self- 
interaction), and v xc (r) is a local exchange-correlation 
potential, where the self-interaction energy is partially 
cancelled. How do the terms in Eq. ( ||) contribute into 
self-trapping? One of the main factors facilitating self- 
trapping is lattice polarization. In alkali halide crystals, 
the second (and perhaps dominant) negative term is the 
energy gain due to the formation of the cr-bond (localised 
doubly occupied orbital) in the molecular ion (X = 
F, CI, Br, I). The factors favouring self-trapping are par- 
tially counterbalanced by an increase in the kinetic en- 
ergy of localizing electrons. The latter results in a split- 
ting of corresponding single particle levels from the edge 
of the valence band. First, we note that any perturba- 
tion in the external potential, v(r), alone cannot cause 
self-trapping, since it does not depend on the number of 
electrons or their state, and self-trapping does not occur 
in a neutral crystal in its ground state. This is in con- 
trast to polaron localization by defects, where the v(r) 
contribution to localization is dominant. Therefore the 
critical localization terms for self-trapping are contained 
in the Hartree (electronic polarisation) and exchange- 
correlation (bond formation) terms. 

Apparently, an accurate calculation of these contribu- 
tions presents a serious problem for the KS DFT, for the 
following reasons: 

1) First, let us consider the formation of the quasi- 
molecule. For that purpose we calculate the dissociation 
curve for the free Cl^ molecular ion using different den- 
sity functionals and the Gaussian-98 packaged. We shall 
use as a reference the energy curve obtained by the cou- 
pled cluster method with double substitutions (CCD). 
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FIG. 6: The dissociation curve for the Cl^ molecular ion 
calculated using various density functionals (Edis = E(Cl~) + 
E cl o ) . All the local basip calculations are done using the G86- 
311 Gaussian basis setHil. In the plane wave GGA calculations 
the energy cut-off 285.3 eV is used. 



In Fig we observe that both the Unrestricted Hartree- 
Fock and all studied DFT functionals do not predict thed 
correct dissociation behaviour for this open shell ion. 
Moreover, local as well as hybrid density functionals (but 
not the Hartree-Fock method) violate the energy varia- 
tional principle in that the predicted dissociation curves 
asymptotically tend to energies significantly lower than 
the corresponding dissociation limit (i.e. the sum of the 
ground state energies of the isolated negative chlorine 
ion and the neutral chlorine atom). The implication of 
this drawback for the self- trapping is clear - the energy 
gain associated with a decrease of the inter-ionic distance 
in the Cl^ from ~ 4 A (as in a perfect lattice) to the 
T4-center equilibrium configuration (*~ 3 A), is substan- 
tially underestimated by the DFT and overestimated by 
the UHF method. 

Incorrect dissociation behaviour of radical ions in DFT 
was first reported by Bally and SastrjH. The subse- 
quent comprehensive-Analysis of the molecular ion 
by Zhang and Yangeil revealed that the self-interaction 
error in local density functional for such systems becomes 
uncontrollably large. It can be demonstrated that the 
self-interaction error, associated with the single particle 
states, is positive and increases with their localization. 
Therefore, self-interaction, being minimized during self- 
consistency, will also lead to the less localized single par- 
ticle states. 

2) The polarization properties of a many electron sys- 
tem in the density functional theory are related to the 
exchange-correlation kernel^ K xc (r, r') = ^^^f^ry- 
The exact exchange correlation kernel is shown to possess 
the so-called ultra-non-locality property, i.e. the diago- 



nal element of its Fourier transform K(q, q) possesses a 
0(l/q 2 ) divergence at small swavenumbers q. However, as 
first demonstrated in refsH3'E2l, both LDA and GGA lead 
to K xc (rr') which does not have the required 0(1/ q 2 ) 
divergence. One of the implications of the non-divergent 
exchange-correlation kernel in KS DFT is that the macro- 
scopic polarizability of the system (calculated as a KS 
polarizability) depends on the number of electrons and 
on the occupancy of the KS bands. That is, a calculated 
KS polarizability of the dielectric will be generally differ- 
ent in a singlet and triplet state, and for a system with or 
without the hole (both systems are technically metallic). 
Thus, the reliability of polarization contribution in KS 
DFT for the triplet or charged systems is also essentially 
uncontrolled. 

3) As mentioned earlier, the largely single particle 
character of the self-trapping problem implies that the 
total spin density may be represented by a single particle 
density of just one orbital. It is known however, that KS 
orbitals do not necessarily carry this physical meaning. 
In particular, by virtue of construction, the KS formalism 
yields the single particle eigenstates which minimize thp 
kinetic energy (as calculated on single particle orbitals )li3. 
Thus, the kinetic energy contribution is always underes- 
timated, leading generally to the least localized set of 
orbitals for a given density. Another delocalizing factor 
is the self-interaction error in the Hartree term, which 
is only partially cancelled in the local approximations to 
the exchange-correlation potentials. Apparently, a pos- 
itive self-interaction grows uncontrollably as the state 
becomes more and more localized. Therefore, any self- 
consistent procedure minimizing the total energy, also 
minimizes self-interaction by delocalizing the orbitals. It 
should be noted, that the effects of minimization of the 
kinetic energy and self-interaction in KS method have 
the same effect and cannot be separated, since any form 
of non-local self-interaction correction leads to the single 
particle equations distinct from those of KS. 

Conceivably the depth of the self-trapping potential in 
polar dielectrics is unlikely to exceed 1 eV. It seems, that 
the spurious delocalizing contributions inherent in the KS 
method always outgrow the localizing potential, so the 
self-trapped state cannot be predicted in this approach. 

This problem is partially resolved by use of hybrid den- 
sity functionals. We have demonstrated for the holes that 
the self-trapped solution gradually becomes energetically 
more stable as the amount of the exact exchange is in- 
creased. However, the conventional hybrid functionals, 
such as B3LYP, failed to predict a stable V^-center, let 
alone its localization energy. This raises the question of 
parametrisation strategies for semiempirical functionals, 
which are not clear in the context of self-trapping. 

Our calculations also highlighted another important 
methodological issue related to the calculations of the 
PES for polar systems within the periodic boundary con- 
ditions. The artificial multipole electrostatic interactions 
between supercells may drastically affect the shape of the 
PES, if the leading multipoles of the charge density vary 
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along the chosen adiabatic coordinate. The dipole-dipole 
interaction in the triplet exciton modelled in a small su- 
percell resulted in the appearance of the energy minimum 
for the STE configuration reported by Perebeinos et alu. 
Interestingly, the authors of refB reported that no stable 
configuration could be found for the Vfc-center. These re- 
sults can be understood by recalling that in the Vfc-center 
calculations the leading term in the electrostatic interac- 
tion between the supercells is the monopole-monopole 
one. This term does not vary with the adiabatic coordi- 
nate and therefore does not affect the shape of the PES. 

Although artificial derealization in KS based DFT 
theory can be expected, it is difficult to assess a priori 
when a specific functional will fail in practical calcula- 
tions. We suggest, that the problem with the localized 
states in KS method will occur in all cases where the 
localization potential is provided mainly by the electron- 
electron term. These include self-trapping phenomena, 
Jahn- Teller systems, low dimensional systems (e.g. Pierls 
instability), and systems with weak impurity perturba- 
tions. A list (by no means complete) of problematic 
DFT calculations where an artificial derealization was 
reported includes: the hole localization in CaO doped 
with Li andpWa and the F-center in LiF as reported by 
Dovesi et alS^b; dimerization in C4jy+a cacban ringsS-l; the 
structure of the Al defect in a-quartzaOl^j-exciton self- 
trappping in presence of thermal disordeiOO; the phase 
diagram of crystallinc-JjUi.-(jrelated to a localization of /- 
type atomic orbitals)E3HEj. 



On the other hand, the use of the Hartree-Fock based 
approaches is also limited since-jaeglecting electron cor- 
relation favours localized statesE9. One should also note 
that frequently used approaches based on tight binding 
methods do not fully account for the electron kinetic en- 
ergy, and thus are also intrinsically biased towards local- 
ized solutions. As seen from these arguments, the DFT 
approach in KS formalism is unable to reliably resolve 
between localized and delocalized situations in systems 
where both type of states are possible. A more adequate 
many-electron theory for the self-trapping problem must 
go beyond the KS methods. Different forms of the den- 
sity functionals free of the self-interaction error should 
be examined. 
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